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We build a sequence of empirical measures on the space B(R+ , R d ) of Revalued cadlag functions 
on R + in order to approximate the law of a stationary Revalued Markov and Feller process (X t ). 
We obtain some general results on the convergence of this sequence. We then apply them to 
Brownian diffusions and solutions to Levy-driven SDE's under some Lyapunov-type stability 
assumptions. As a numerical application of this work, we show that this procedure provides an 
efficient means of option pricing in stochastic volatility models. 
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1. Introduction 

1.1. Objectives and motivations 

In this paper, we deal with an Revalued Feller Markov process (X t ) with semigroup 
(Pt)t>o and assume that (X t ) admits an invariant distribution vq. The aim of this work is 
to propose a way to approximate the whole stationary distribution ¥ Ua of (X t ). More pre- 
cisely, we want to construct a sequence of weighted occupation measures (j/ n ) (w, da)) n >i 
on the Skorokhod space D(R+,R d ) such that n Z^ a j F(a)F Vo (da) a.s. for a 

class of functionals F : D(R + , M d ) which includes bounded continuous functionals for the 
Skorokhod topology. 

One of our motivations is to develop a new numerical method for option pricing in sta- 
tionary stochastic volatility models which are slight modifications of the classical stochas- 
tic volatility models, where we suppose that the volatility evolves under its stationary 
regime. 



This is an electronic reprint of the original article published by the ISI/BS in Bernoulli, 
2009, Vol. 15, No. 1, 146-177. This reprint differs from the original in pagination and 
typographic detail. 



1350-7265 © 2009 ISI/BS 



Approximation of the distribution of a stationary Markov process 



147 



1.2. Background and construction of the procedure 

This work follows on from a series of recent papers due to Lamberton and Pages ([12, 13]), 
Lemaire ([14, 15]) and Panloup ([18, 19, 20]), where the problem of the approximation 
of the invariant distribution is investigated for Brownian diffusions and for Levy-driven 
SDE's. 1 In these papers, the algorithm is based on an adapted Euler scheme with de- 
creasing step (7fe)fe>i. To be precise, let (F„) be the sequence of discretization times: 
T = 0, r n = J^ILi 7fc for every n > 1, and assume that T n — ► +oo when n — > +oo. Let 
(Xr rl ) n >o be the Euler scheme obtained by "freezing" the coefficients between the T„'s 
and let (//„)„>! be a sequence of positive weights such that H n := ^fc=i Vk ~ * +°o when 
k — > +oo. Then, under some Lyapunov-type stability assumptions adapted to the stochas- 
tic processes of interest, one shows that for a large class of steps and weights (?7„,7 n ) n >i, 

1 " - — f 

v n {ujJ)-.= —^r) k f{X Tk _ 1 ) n ^° I f{x)v n {&x) a.s., (1) 

Un k=i J 

(at least) 2 for every bounded continuous function /. 

Since the problem of the approximation of the invariant distribution has been deeply 
studied for a wide class of Markov processes (Brownian diffusions and Levy-driven SDE's) 
and since the proof of (1) can be adapted to other classes of Markov processes under some 
specific Lyapunov assumptions, we choose in this paper to consider a general Markov pro- 
cess and to assume the existence of a time discretization scheme (^r fc )fc>o such that (1) 
holds for the class of bounded continuous functions. The aim of this paper is then to inves- 
tigate the convergence properties of a functional version of the sequence (p n (uj ,&a)) n >\. 

Let (X t ) be a Markov and Feller process and let (X t )t>o be a stepwise constant time 
discretization scheme of (X t ) with non-increasing step sequence (7 n )n>i satisfying 

n 

lim 7„ = 0, T„ := V" j k "n±>°° +oc . (2) 

// hoc * — * 

k=l 

Letting To := and Xq = xq £ M. d , we assume that 

X t =Xr n Vie [L„,L, 1+ i[ (3) 

and that (Xr n ) n >o can be simulated recursively. 

We denote by (Tt)t>a and {Ft)t>Q the usual augmentations of the natural filtrations 
(a(X s ,0 < s < t)) t >o and (a(X s ,0 < s < t)) t >o, respectively. 

1 Note that computing the invariant distribution is equivalent to computing the marginal laws of the 
stationary process (Xt) since v^Pt = vq for every t > 0. 

2 The class of functions for which (1) holds depends on the stability of the dynamical system. In 
particular, in the Brownian diffusion case, the convergence may hold for continuous functions with 
subexponential growth, whereas the class of functions strongly depends on the moments of the Levy 
process when the stochastic process is a Levy-driven SDE. 
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For k > 0, we denote by {X[ k ^)t>o the shifted process defined by 

A t -- A r fc +t- 

In particular, X^ = X t . We define a sequence of random probabilities {u^ n \u!,da)) n >i 
on P(M + ,R d ) by 



1 n 

v {n) (w, da) = — Vkl { x(k-D 



" k=l 

where {rjk)k>i is a sequence of weights. For t > 0, (i/j (w, dx))„>i will denote the se- 
quence of "marginal" empirical measures on M. d defined by 

1 ™ 

k—1 

1.3. Simulation of (o^F))^ 

For every functional F : B(R+, R d ) -> R, the following recurrence relation holds for every 
n>l: 

^ n+1) {cu,F) = u< n) (u>,F) + %^-{F{X^{u)) - v ^\oj,F)). (4) 

Then, if T is a positive number and F : B(R + , R d ) — > R is a functional depending only 
on the trajectory between and T, {v^ n > {u , F)) n >\ can be simulated by the following 
procedure. 

Step 0. (i) Simulate {X^° )t>o ° n [O^L that is, simulate {Xr k )k>o for ft = 
0,...,N{0,T), where 

JV(n, T) := inf{ft > n, T k+1 - T n > T} 

= max{fe>0,r fc -r„<T}, n>0,T>0. 

Note that n t— > N(n,t) is an increasing sequence since (7„) is non-increasing, and that 

^N( n ,T) — r„ < t < r N ( n ,t)+i — r„. 

(ii) Compute F((X t (0) ) t > ) and i/W^.F). Store the values of (X T „) for ft = 
1,...,JV(0,T). 

Step n (n > 1). (i) Since the values {X-p k )k>o are stored for ft = n, . . . , iV(n — 1, T), 
simulate {Xr k )k>o for ft = N(n— 1, T) + 1, . . . , N{n, T) in order to obtain a path of {XI ) 
on [0,T]. 

(ii) Compute F((X t (n) ) t>o) and use (4) to compute z^" +1 ' F). Store the values of 
(X r J for fc = n + l,...,#(n,T). 
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Remark 1. As shown in the description of the procedure, one generally has to store 
the vector [-Xr„, ■ ■ ■ , X^ N{n T) ] at time n. Since (7„) is a sequence with infinite sum that 
decreases to 0, it follows that the size of this vector increases "slowly" to +oo. For 
instance, if 7„ = Cn~ p with p £ (0, 1), its size is of order n p . However, it is important 
to remark that even though the number of values to be stored tends to +oo, that is 
not always the case for the number of operations at each step. Indeed, since X( n+1 *> 
is obtained by shifting X^ n \ it is usually possible to use, at step n+1, the preceding 
computations and to simulate the sequence (i ;l (A > ( n )))„>o in a "quasi-recursive" way. 
For instance, such remark holds for Asian options because the associated pay-off can be 
expressed as a function of an additive functional (see Section 5 for simulations). 

Before outlining the sequel of the paper, we list some notation linked to the spaces 
B(R + ,R d ) and B([0,T],R d ) of cadlag Revalued functions on R + and [0,T], respectively, 
endowed with the Skorokhod topology. First, we denote by d\ the Skorokhod distance 
on B([0, l],R d ) defined for every a, [3 £ ©([0, l],R d ) by 

d\ («,/?)= inf < maxl sup \a(t) — [3(X(t))\, sup 
AeAl I \te[o,i] o<s<t<i 

where Ai denotes the set of increasing homeomorphisms of [0,1]. Second, for T > 0, 
4> T :B(R + ,R d ) i-> ©([0, l],R d ) is the function defined by (</> T (a))(s) = a(sT) for every s £ 
[0, 1]. Wc then denote by d the distance on B(R+,R d ) defined for every a, (3 £ B(R+,R d ) 

by 

r+oo 

d(a,f3)= c- t (lAd 1 (0 t (a),0 t (/3)))dt. (6) 

JO 

We recall that (D(R+,R d ), d) is a Polish space and that the induced topology is the usual 
Skorokhod topology on D(R + ,R d ) (see, e.g., Pages [16]). For every T > 0, we set 

T> T = f| <7(7r u ,0<M<s), 

s>T 

where ir s : D(R+, R d ) -> R d is defined by ir s {a) = a(s). For a functional F:B(R + ,R d ) -> 
R, Ft denotes the functional defined for every a £ D(R+,R d ) by 

F T {a) = F(a T ) with a T (t) = a(t A T) Vi > 0. (7) 

Finally, we will say that a functional -F:B(R+,R d ) — > R is S'fc-continuous if F is contin- 

(Sk) 

uous for the Skorokhod topology on B(R + ,R rf ) and the notation "=>" will denote the 
weak convergence on D(R + ,R rf ). 

In Section 2, we state our main results for a general Revalued Feller Markov process. 
Then, in Section 3, we apply them to Brownian diffusions and Levy-driven SDE's. Section 
4 is devoted to the proofs of the main general results. Finally, in Section 5, we complete 
this paper with an application to option pricing in stationary stochastic volatility models. 
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2. General results 

In this section, we state the results on convergence of the sequence (p( n >(u>, da)) n >i when 
(X t ) is a general Feller Markov process. 

2.1. Weak convergence to the stationary regime 

As explained in the Introduction, since the a.s. convergence of (pq (lu, dx))„>i to the 
invariant distribution vq has already been deeply studied for a large class of Markov 
processes (Brownian diffusions and Levy driven SDE's), our approach will be to derive 
the convergence of (i/ n ) (to, da))„>i toward P„ from that of (i/q (to, dx))„>i to the 
invariant distribution vq. More precisely, we will assume in Theorem 1 that 
(Co.i): (Xt) admits a unique invariant distribution Vq and 

Vq {u!,ax) =£- ^o(da;) a.s., 

whereas in Theorem 2, we will only assume that 
(Co, 2): {vq \ijj,dx)) n >i is a.s. tight on R d . 

We also introduce three other assumptions, (Ci), (C 2 ) and (C 3j£ ), regarding the conti- 
nuity in probability of the flow x 1— ► (X®), the asymptotic convergence of the shifted time 
discretization scheme to the true process (X t ) and the steps and weights, respectively. 

(Ci): For every a; G M d , e > and T > 0, 

lim sup P ( sup \Xf - X't a I > e ) = 0. (8) 

xq^x \0<t<T J 

(C2): (Xt) is a non-homogeneous Markov process and for every n > 0, it is possible to 
construct a family of stochastic processes (Y t ) x gRcj such that 

(i) C(YM) D(R ^' Rd) C(X^\X { n) = x); 

(ii) for every compact set K of R d , for every T > 0, 

sup sup |y t (n * a) -jrf| n =^ oo in probability. (9) 

x£K0<t<T 

(C 3 , e ): For every n > 1, ry„ < C^ n R e n . 

Remark 2. Assumption (C2) implies, in particular, that asymptotically and uniformly 
on compact sets of M. d , the law of the approximate process (X^), given its initial value, 
is close to that of the true process. 

If there exists a unique invariant distribution vq, the second part of (C2) can be relaxed 
to the following, less stringent, assertion: for all e > 0, there exists a compact set A e C R d 
such that iso(A^) < e and such that 

sup sup \Y t (n ' x) - Xf I n =^°° in probability. (10) 

igA, 0<t<T 
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This weaker assumption can some times be needed in stochastic volatility models like 
the Heston model (see Section 5 for details). 

The preceding assumptions are all that we require for the convergence of (y( n ' [to, da))„>i 
to P„ along the bounded S'fc-continuous functionals, that is, for the a.s. weak conver- 
gence on B(R + ,R d ). However, the integration of non-bounded continuous functionals 
F : B([0, T], M. d ) — > R will need some additional assumptions, depending on the stability 
of the time discretization scheme and on the steps and weights sequences. We will sup- 
pose that F is dominated (in a sense to be specified later) by a function V : M. d — > R + 
that satisfies the following assumptions for some s > 2 and e < 1. 

H(s,e): For every T > 0, 



(i) sup I 



n>l 



sup V s (Y t {n ' x) ) 



0<t<T 



<C T V s (x), 



(ii) sup^" ) (V)<+oo, 



(hi) ^^E[V 2 (Xr fc _ 1 )]<+oo, 
fc>i k 

(iy) J^^^-nV^HXr^K+oc, 

k>l k 

where T i-> C T is locally bounded on M.+ and AN(k, T) = N(k, T) - N(k - 1, T). 
For every e < 1, we then set 

JC{e) = {Ve C(M' i ,R + ),H(s,e) holds for some s > 2}. 



Remark 3. Apart from assumption (i), which is a classical condition on the finite time 
horizon control, the assumptions in H(s,e) strongly rely on the stability of the time 
discretization scheme (and then, to that of the true process). More precisely, we will sec 
when we apply our general results to SDE's that these properties are some consequences 
of the Lyapunov assumptions needed for the tightness of {v\^ {uj,dx)) n >i. 

We can now state our first main result. 

Theorem 1. Assume (Co.i), (Ci), (C2) and (C3. e ) with eg (— 00, 1). Then, a.s., for 
every bounded Sk-continuous functional i^:B(R + ,R d ) — > R, 

"=±>°° J F(a)¥ V0 (da), (11) 

where P„ denotes the stationary distribution of (X t ) (with initial law Uq). 

Furthermore, for every T > 0, for every non-bounded Sk-continuous functional 
J F:B(R+,K d ) -► R, (11) holds a.s. for F T (defined by (7)) if there exists V G JC{e) and 
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p G [0, 1) such that 

\F T (a)\<C sup V p {a t ) Va G D(R+, R d ). (12) 

0<t<T 

In the second result, the uniqueness of the invariant distribution is not required and 
the sequence (vq {u>,dx)) n >i is only supposed to be tight. 

Theorem 2. Assume (Co^), (Ci), (C2) and (C3 i£ ) with £G (— 00, 1). Assume that 
(vq (w,da;))„>i is a.s. tight on . We then have the following. 

(i) The sequence (v( n > (cj, da)) n >i is a.s. tight on U>(R + ,R d ) and a.s., for ev- 
ery convergent subsequence (n/ c (w))„>i ; for every bounded Sk-continuous functional 
F:B(R + ,R d ) —>R, 

f F{a)W Voo (da), (13) 



where Pj,^ is the law of {Xt) with initial law being a weak limits for {uq {u>,dx)) n >i. 

Furthermore, for every T > 0, for every non-bounded Sk-continuous functional 
F:B(R+,R d ) -> R, (13) holds a.s. for F T if (12) is satisfied with V G fC{e) and p G [0, 1). 
(ii) If, moreover, 

iV«^"^0, (14) 

then Voo is necessarily an invariant distribution for the Markov process {Xt). 

Remark 4- Condition (14) holds for a large class of steps and weights. For instance, 
if T] n = C\n~ pi and -f n = C-2n~ p ' 2 with pi G [0, 1] and P2 G (0, 1], then (14) is satisfied if 
pi = or if pi G (max(0, 2p 2 — 1), 1). 



2.2. Extension to the non-stationary case 

Even though the main interest of this algorithm is the weak approximation of the pro- 
cess when stationary, we observe that when v is known, the algorithm can be used to 
approximate P^ if po is a probability on R d that is absolutely continuous with respect 

to Vq. 

Indeed, assume that pg{dx) = (f){x)vQ{dx) , where 0:R d ^R is a continuous non- 
negative function. For a functional F : D(R + , R d ) — » R, denote by F^ the functional de- 
fined on D(R+,R d ) by F^{a) = F(a)0(a(O)). 

Then, if v( n \u>,da) P„ (da;) a.s., we also have the following convergence: a.s., for 
every bounded S k- continuous functional F:D(R + ,R d ) — > R, 

v W(w,F+) n ^?° J F+(a)W V0 (da)= J F{a)P^ {da). 
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3. Application to Brownian diffusions and 
Levy-driven SDE's 

Let (X t ) t >o be a cadlag stochastic process solution to the SDE 

dX t = b(X t -)dt + a(X t -)dW t + K(X t -)dZ t , (15) 

where b : R d -» R d , a : R d i-» M d . e (set of d x I real matrices) and n : R d i-» M d;£ are con- 
tinuous functions with sublinear growth, (Wt)t>o is an ^-dimensional Brownian motion 
and (Z t )t>o is an intcgrable purely discontinuous R^-valued Levy process independent of 
(Wt)t>o with Levy measure n and characteristic function given for every t > by 

E[ e i < u '^>]=exp 

Let ( / y n )n>i be a non-increasing step sequence satisfying (2). Let {U n ) n >\ be a sequence 
of i.i.d. random variables such that U± =Af(0,Ie) and let £ := (£ n ) n >i be a sequence of 
independent R^-valued random variables, independent of (U n ) n >x- We then denote by 
(X t )t>o the stepwise constant Euler scheme of (X t ) for which (Xr„) n >o is recursively 
defined by X = x G M d and 

X r „ +1 =X Fn + 7)l+ ife(X r J + V7^Tcr(X r J[/ n+ i + K (X r J£„+i. (16) 

We recall that the increments of (Z t ) cannot be simulated in general. That is why we 
generally need to construct the sequence (£„) with some approximations of the true 
increments. We will come back to this construction in Section 3.2. 

As in the general case, we denote by (X^) k > and {v^ (u, da)) n >i the sequences of 
associated shifted Euler schemes and empirical measures, respectively. 

Let us now introduce some Lyapunov assumptions for the SDE. Let £Q(R d ) denote 
the set of essentially quadratic C 2 -functions V :R d — > R+ such that limV(x) = +oo as 
\x\ — > +oo, | VV| < CVV and D 2 V is bounded. Let a £ (0, 1] denote the mean reversion 
intensity. The Lyapunov (or mean reversion) assumption is the following. 

(S a ): There exists a function V £ £Q(R d ) such that: 

(i) \b\ 2 <CV a , Ti(aa*(x)) + \\K(x)\\ 2 W = + °°o(r(i)); 

(ii) there exist [3 G R and p > such that (W, b}</3- pV a . 

From now on, we separate the Brownian diffusions and Levy-driven SDE cases. 
3.1. Application to Brownian diffusions 

In this part, we assume that k = 0. We recall a result by Lambcrton and Pages [13]. 

Proposition 1. Let a € (0, 1] such that (S a ) holds. Assume that the sequence (7y n /7n)n>i 
is non-increasing. 



t( / e i{u ' v) - l-i(u,y)ir(dy) 
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(a) Let (9 n )n>i be a sequence of positive numbers such that J2n>i® n "f n < +°° an ^ 
that there exists no G N such that (9 n )n>n * s non-increasing. Then, for every positive r, 



(b) For every r > 0, 



^e n7n E[y r (Xr„_ 1 )]<+oo. 



+00 a.s. (17) 

71>1 



Hence, the sequence (y^ (u),dx)) n >i is a.s. tight. 

(c) Moreover, every weak limit of this sequence is an invariant probability for the SDE 
(15). In particular, if (Xt)t>o admits a unique invariant probability vq, then for every 
continuous function f such that f < CV r with r > 0, linin^oo z/q (lo, f) = vo{f) a.s. 

Remark 5. For instance, if V(x) = 1 + |x| 2 , then the preceding convergence holds for 
every continuous function with polynomial growth. According to Theorem 3.2 in Lemairc 
[14] , it is possible to extend these results to continuous functions with exponential growth, 
but it then strongly depends on a. Further the conditions on steps and weights can be 
less restrictive and may contain the case rj n = 1 , for instance (see Remark 4 of Lamberton 
and Pages [13] and Lemaire [14]). 

We then derive the following result from the preceding proposition and from Theorems 
1 and 2. 

Theorem 3. Assume that b and a are locally Lipschitz functions and that k = 0. Let 
a £ (0, 1] such that (S a ) holds and assume that (r] n /j n ) is non-increasing. 

(a) The sequence (c'"'(w,da)) n >i is a.s. tight on C(R+,R d ) 3 and every weak limit 
of (i»(w,da)) „>i is the distribution of a stationary process solution to (15). In par- 
ticular, when uniqueness holds for the invariant distribution Vq, a.s., for every bounded 
continuous functional F :C(M.+ ,W i ) — > R, 



J F(x)W„ (dx). (18) 

(b) Furthermore, if there exists s £ (2, +00) and no £ N such that 

f AN(k,T) \ jSr AN(k,T) ^ 
zs non-increasing and > < +00, (19) 

V ikH s k y„> no ^ HI 



3 C(E + ,R d ) denotes the space of continuous functions on R_|_ with values in R d endowed with the 
topology of uniform convergence on compact sets. 
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then, for every T > 0, for every non-bounded continuous functional F:C(M. + ,R d ) — ► R, 

(18) holds for Ft if the following condition is satisfied: 

3r>0 suchthat \F T (a)\ <C sup V r (a t ) Va e C(M + , M d ). 

Q<t<T 

Remark 6. If rj n = C\n~ Pl and j n = C-2n~ P2 with < p2 < pi < 1, then for s <G (1, +oo), 

(19) is fulfilled if and only if s > 1/(1 — p\). It follows that there exists s e (2,+oo) such 
that (19) holds as soon as p\ < 1. 

Proof of Theorem 3. We want to apply Theorem 2. First, by Proposition 1, assumption 
(Co, 2) is fulfilled and every weak limit of (v^ (w, da;)) is an invariant distribution. Second, 
it is well known that (Ci) and (C2) are fulfilled when b and a are locally Lispchitz 
sublinear functions. Then, since (C3 iB ) holds with e = 0, (18) holds for every bounded 
continuous functional F. Finally, one checks that H(s,0) holds with V :— V r (r > 0). 
It is classical that assumption (a) is true when b and a arc sublinear. Assumption (b) 
follows from Proposition 1(b). Let 9 n ,i = ?7 n /(7„i^) and 9 n .2 = AiV(n,T)/(7 n H*). Using 
(19) and the fact that (rj n /'-f n ) is non-increasing yields that (# n ,i) and {0 n ^) satisfy the 
conditions of Proposition 1 (see (35) for details). Then, (iii) and (iv) of H(s,0) arc 
consequences of Proposition 1(a). This completes the proof. □ 

3.2. Application to Levy-driven SDE's 

When we want to extend the results obtained for Brownian SDE's to Levy-driven SDE's, 
one of the main difficulties comes from the moments of the jump component (see Panloup 
[18] for details). For simplification, wc assume here that (Z t ) has a moment of order 
2p > 2, that is, that its Levy measure n satisfies the following assumption with p > 1: 

(Hi): [ n(dy)\y\^<+^. 

We also introduce an assumption about the behavior of the moments of the Levy measure 
at 0: 

(H=): / 7r(d»)|y| 2 «<+oo ) ge[0,l]. 
J\y\<i 

This assumption ensures that (Z t ) has finite 2g-variations. Since f\ y \<x \y\ 27r (dy) is finite, 
this is always satisfied for q=l. 

Let us now specify the law of (£„) introduced in (16). When the increments of (Zt) can 
be exactly simulated, we denote by (E) the Euler scheme and by (£u,e) the associated 
sequence 

£,n,E = Z ln Vn>l. 
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When the increments of (Z t ) cannot be simulated, we introduce some approximated Eulcr 
schemes (P) and (W) built with some sequences (£n,p) and (£ n ,w) 01 approximations of 
the true increment (see Panloup [19] for more detailed presentations of these schemes). 
In scheme (P), 



where (Z. in ) n >i a sequence of compensated compound Poisson processes obtained by 
truncating the small jumps of (Z t )t>o- 



where (u n ) n >i is a sequence of positive numbers such that u n — ► 0. We recall that 
Z. n n Zl^° z locally uniformly in L? (see, e.g., Protter [21]). 

As shown in Panloup [19], the error induced by this approximation is very large when 
the local behavior of the small jumps component is irregular. However, it is possible to 
refine this approximation by a Wienerization of the small jumps, that is, by replacing 
the small jumps by a linear transform of a Brownian motion instead of discarding them 
(see Asmussen and Rosinski [2]). The corresponding scheme is denoted by (W) with £, n ,w 
satisfying 

in,W = €n,P + y/lZQnh-n Vn > 1, 

where (A ra )„>i is a sequence of i.i.d. random variables, independent of (£ n ,p)n>i and 
(U n ) n >i, such that Ai =Jx(Q,Ii) and (Q n ) is a sequence of i x I matrices such that 



We recall the following result obtained in Panloup [18] in our slightly simplified frame- 
work. 

Proposition 2. Let a G (0, 1], p > 1 and q £ [0, 1] such that (Hi), (H^) and (S a ) hold. 

Assume that the sequence {r]n/ln)n>i is non- increasing. Then, the following assertions 
hold for schemes (E), (P) and (W). 

(a) Let {9 n ) satisfy the conditions of Proposition 1. Then, ^2 n>1 6n'JnM[V p+a ~ 1 (Xr„- 1 )] < 





(20) 




+oo. 



(b) We have 



n>l 



( n \ UJ ,V^ 2+a - 1 )<+. 



DO 



a.s. 



(21) 



Hence, the sequence {v^ (u>,dx)) n >i is a.s. tight as soon as p/2 + a — 1 > 0. 
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(c) Moreover, if Ti(cra*) + \\k\\ 2 ^ < CV p ^ 2+a - 1 , then every weak limit of this sequence 
is an invariant probability for the SDE (15). In particular, if (X t )t>o admits a unique 
invariant probability vq, for every continuous function f such that f = o{V p / 2+a ~ 1 ) , 
Ym\ n ^ 00 v^ l \uj,f) = h>u{f) a.s. 

Remark 7. For schemes (E) and (P), the above proposition is a direct consequence of 
Theorem 2 and Proposition 2 of Panloup [18]. As concerns scheme (W), a straightforward 
adaptation of the proof yields the result. 



Our main functional result for Levy-driven SDE's is then the following. 



Theorem 4. Let a € (0, 1] and p>l such that p/2 + a — 1 > and let q € [0, 1] . Assume 
(Hp), (Hq) and (S a ). Assume that b, a and k are locally Lipschitz functions. If, more- 
over, (?7n/7n)n>i is non-increasing , then the following result holds for schemes (E), (P) 
and (W). 

(a) The sequence (y( n > (lj, da)) ra >i is a.s. tight on D(K + ,l <i ). Moreover, if 

Tr( £ 7 t 7') + ||K|| 2 9<Cy/ 2+ °- 1 or — V max ^ 0, (22) 

Ha f-('>fc+l 7<-l 

K — 1 

then every weak limit of (i/ n ) (u>, da))„>i is the distribution of a stationary process solu- 
tion to ( 15). 

(b) Assume that the invariant distribution is unique. Let e < such that (C3. e ) holds. 
Then, a.s., for every T > 0, for every Sk-continuous functional F : D(R+, R d ) -> E, (18) 
holds for Ft if there exist p £ [0,1) and s > 2, such that 

\F T (a)\<C sup V {p{p+a - 1))/s (a t ) Vae©(R+,R d ) 

0<t<T 

and if 

fAN(k,T)\ 1 ^AN(k,T) 

I s(l-e) J lS non ~ mcreasm 9 an d 2_^i g(l-e) ^ +°°- (23) 

\lkH S k 6 ) n>n fe >! Hi 



Remark 8. In (22), both assumptions imply the invariance of every weak limit of 
(vq (w,dx)). These two assumptions arc very different. The first is needed in Proposition 
2 for using the Echeverria- Weiss invariance criteria (see Ethier and Kurtz [7], page 238, 
Lamberton and Pages [12] and Lemaire [14]), whereas the second appears in Theorem 
2, where our functional approach shows that under some mild additional conditions on 
steps and weights, every weak limit is always invariant. 

For (23), we refer to Remark 6 for simple sufficient conditions when (7„) and (q n ) are 
some polynomial steps and weights. 
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4. Proofs of Theorems 1 and 2 

We begin the proof with some technical lemmas. In Lemma 1, we show that the a.s 
weak convergence of the random measures ((/'"'(w,da)) n >i can be characterized by the 
convergence (11) along the set of bounded Lipschitz functionals F for the distance d. 
Then, in Lemma 2, we show with some martingale arguments that if the functional 
F depends only on the restriction of the trajectory to [0,T], then the convergence of 
(i/ n ) (u>, F)) n >\ is equivalent to that of a more regular sequence. This step is fundamental 
for the sequel of the proof. 

Finally, Lemma 4 is needed for the proof of Theorem 2. We show that under some mild 
conditions on the step and weight sequences, any Markovian weak limit of the sequence 
(y'"'(w,da)) n >i is stationary. 

4.1. Preliminary lemmas 

Lemma 1. Let (E,d) be a Polish space and let V(E) denote the set of probability 
measures on the Borel a -field B(E), endowed with the weak convergence topology. Let 
(/i'"'(w,da))„>i be a sequence of random probabilities defined on Q x B(E). 

(a) Assume that there exists /i*- 00 ) G V(E) such that for every bounded Lipschitz func- 
tion F:E^R, 

nW{u,F) n ^±?° /iW(f) a.s. (24) 

Then, a.s., (/i'"'(w,da)) n >i converges weakly to //°°) onV(E). 

(b) Let hi be a subset ofV(E). Assume that for every sequence (i 7 fc)fc>i of Lipschitz 
and bounded functions, a.s., for every subsequence (n^"^ n '' (lj, da)), there exists a sub- 
sequence (/i^"°^"("))(a;,da)) and a U -valued random probability (j,(°°> (uj, da) such that 
for every k>l, 

^' Mn) \u J ,F k ) n ^t a ^\uj,F k ) a.s. (25) 
Then, (fx^ n ' (w, da)) n >i is a.s. tight with weak limits inlA. 

Proof. We do not give a detailed proof of the next lemma, which is essentially based 
on the fact that in a separable metric space (E,d), one can build a sequence of bounded 
Lipschitz functions (g k )k>i such that for any sequence (fJ. n )n>i of probability measures 
on B(E), (/i n )n>i weakly converges to a probability fi if and only if the convergence 
holds along the functions g k , k>l (see Parthasarathy [22], Theorem 6.6, page 47 for a 
very similar result). □ 

For every n > 0, for every T > 0, wc introduce r(n,T) defined by 

r(n, T) := min{fc > 0, N(k, T)>n} = min{fc < n, T k + T> T n }. (26) 
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Note that for k € {0, . . . , r(n,T) - 1}, {x[ k \{) <t<T} is ^r n -measurable and 

t - 7 T (n,T)-i < r„ - r T („, T ) < T. 



Lemma 2. Assume (C3 j£ ) with e < 1. Lei F:B(R + ,M d ) — > R be a Sk-continuous func- 
tional. Let {Qk) be a filtration such that Tv k C Qk for every k > 1. Then, for any T > 0: 

(a) z/ -Ft (defined by (7)) is bounded, 

1 " _ 

— ^ ??fe (F T (X( fc - 1 ))-E[F T (X( fe - 1 ))/g fc _ 1 ])"=±> oo o. fl .; (27) 
" fc=i 

(b) if Ft is not bounded, (27) holds if there exists V:R d — >R+, satisfying H(s,e) for 
some s>2, such that \Fr(ct) \ < Csup 0<t<T V(at) for every a E H)(M. + ,M. d ) ; furthermore, 

sup v {n) {u), F T ) < +oo a.s. (28) 
n>l 

Proof. We prove (a) and (b) simultaneously. Let T( fc ) be denned by T^ fe ' = i*V(X( fe '). 
We have 



i^^( T(fe " 1) - E [ T(fe_1) /^-i]) 
Hn fc=i 

= 7rE^( T(fc_1) - E [ T( ^ 1) /^]) (29) 

W " fc=i 
1 " 

+ 7rE' ? ' fe ( E [ T(fe_1) /^] -EIT^" 1 )/^-!]). (30) 

11 n 



k=l 



We have to prove that the right-hand side of (29) and (30) tend to a.s. when n — > +oo. 

We first focus on the right-hand side of (29). From the very definition of r(n,T), we 
have that {X{ , < t < T} is JFr„ -measurable for k G {0, . . . , r(n, T) — 1}. Hence, since 
Ft is cr(7r s ,0 < s < T)-measurable and jFr„ c5 n , it follows that is (J n -measurable 
and that T( fe) = E[T^ /Q n ] for every k < r(n, T) — 1. Then, if Ft is bounded, we derive 
from (Cs. e ) that 



^E^( T( *" 1) - E [ T(fc ~ 1) /^ 



k—l 



™ fc=r(n,T) + l ™ fc=r(n,T) + l 

C 

< !_ r (r n - r r( „ jT) ) 

tin 
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Hn 

where we used the fact that (H n ) n >i and (7 n )n>i are non-decreasing and non-increasing 
sequences, respectively. 

Assume, now, that the assumptions of (b) are fulfilled with V satisfying H(s,e) for 
some s > 2 and e < 1. By the Borel-Cantelli-like argument, it suffices to show that 



E E 



i n 



n k=r(n,T)+l 



< +oo. 



(31) 



Let us prove (31). Let a k := r) { k s ~ 1)/s and b k {u>) := ^"(T^- 1 ) - E[T^-^ /Q n ]). The 
Holder inequality applied with p = s/(s — 1) and <f = s yields 



^ a k b k (uj) 
fe=r(n,r)+l 



< 



E * E •»iT( fc - i) -E[T(*- i )/ft l ]r 

vfc=T(n,T)+l / \fe=r(n,T)+l / 



Now, since Ft (a) < sup 0<t<T V(a), it follows from the Markov property and from 
H(s,e)(i) that 



sup V s {X^)/T Tk 

0<t<T 



E[|F r (xW)|7JT»,]<ci: 

Then, using the two preceding inequalities and (C3 j£ ) yields 

n 

J2 %(T (fc - 1) -E[T( fe - 1 V£„]) 

k=T(n,T) + l 



<C T V s (X Tk ). 



<c[ E * 

\fc=-r(ri,T) + l / 

< c f E ^1 E 



kfc=T(n,T)+l / 



sup V'^r,.,; 

fc=-r(n,T) + l 



\fc=r(ri,T) + l / 



sup v s (x; (n ' T) ) 

,te[0,S(n,T)] 



where S(n,T) = L„_i — T T r n> T) and C does not depend n. By the definition of r(n,T), 
S(n,T) < T. Then, again using H(s,e)(i) yields 



E E 



H' 



fc=r(n,T) 



— ' H" 

„>1 iln 



T))J 
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Since n i— > N(n, T) is an increasing function, n i— > r(n, T) is a non-decreasing function 
and Card{n, r(n,T) = fc} = AJV(fc + l,T) := iV(fe + 1,T) - N(k,T). Then, since n^H n 
increases, a change of variable yields 



71>1 I 



1 



k=r(n,T) + l 



k>l U k 



by H(s,e)(iv). 

Second, we prove that (30) tends to 0. For every n > 1, we let 



M„ = V^(E[T( fe - 1 )/S„]-E[T( fe - 1 )/a fe -i]). 



fe=i 



(32) 



The process {M n ) n >\ is a (C/ n )-martingale and we want to prove that this process is 
^-bounded. Set = E[F T (X^)/g n ] - E[F T (X^)/g k }. Since F T is cr(7r s ,0 < s < 

T) -measurable, the random variable is ^r N(k T) -measurable. Then, for every i g 

{N(k,T), ...,n}, <J> (fc > n) is ^-measurable so that 

E [$(i,n) $ (fc,n)] = E[$(fc.n)]g[$(i.«)/0.]] = O. 

It follows that 

2 N(k-l,T)An 

E[M^]=y 4E[(l (fc - lin) ) 2 ]+2yJ V ^-Et*^- 1 '")*^- 1 '^]. (33) 

fe>l « fe>l K i=fc+l 



Then, 



2 JV(fe-l,T) 

su P E[M2]<V ^supE[(<f>( fc - 1 ^) 2 ]+2y ^ y -^supE^" 1 '™^- 1 ^] 



<^fE^^p E [( $(fe - 1,ri) ) 2 ] 



(34) 



y^L_ y 7i supE[$( < - 1 '")*(*- 1 ' n )] 



2V(fe— 1,T) 



fe>l " *=fc+l 
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where, in the second inequality, we used assumption (C3 j£ ) and the decrease of i 
\jH\~ e . Hence, if Ft is bounded, using the fact that Y^=k+x Ji<T yields 



supE[M2]<C]T-^<C 
fc>i H k 



Hi 



2-e 



du 



'/i 



< +00 



(35) 



since e < 1. Assume, now, that the assumptions of (b) hold and let Ft be dominated 
by a function V satisfying H(s,e). By the Markov property, the Jensen inequality and 
H(s, £ )(i), 



E[($( fe <™)) 2 ] <CE 



E 



sup V 2 (X«)/^ 

0<t<T 



<C T E[V 2 (A r J 



We then derive from the Cauchy-Schwarz inequality that for every n, k > 1 , for every 
i€{k,...,N(k,T)}, 



|E[$(»» $(*,»»)] J < cJE[V 2 (X r ,)}jE[V 2 (Xr k )} < C sup E[V 2 (X^ k) )} < CE[V 2 (X r J], 

v te[o,T] 

where, in the last inequality, we once again used H(s,e)(i). It follows that 
su P E[M 2 ] < CJ2 -J^nV^Xr^)] < +00, 

by H(s,e)(iii). Therefore, (34) is finite and (M ra ) is bounded in L 2 . Finally, we derive 
from the Kronecker lemma that 



1 n 



+ 00 



a.s. 



fe=i 



As a consequence, sup n>1 v^ n > (u) , Ft) < +00 a.s. if and only if 

sup — V E[F T {X^)/T k -i] <+oo a.s. 



fe=l 



This last property is easily derived from H(s,e)(i) and (ii). This completes the proof. □ 

Lemma 3. (a) Assume (Ci) and let xq € R d . We then have lim^^ E[d(X x , X x °)] = 0. 
In particular, for every bounded Lispchitz (w.r.t. the distance d) functional F :B(R+,R rf ) — » 
R, the function $ F defined by $ F (x) = E[F(X X )} is a (bounded) continuous function on 
R d . 

(b) Assume (C2). For every compact set K C R d , 

sup E[d(Y n ' x , X x )] "^i 00 0. (36) 
xeK 
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Set ®%(x) = E[F(Y n > x )} . Then, for every bounded Lispchitz functional F : D(R + , R d ) -> R, 
sup - ™^±?° for every compact set K C R d . (37) 

Proof, (a) By the definition of d, for every a, [3 S B(R+,R d ) and for every T > 0, 

< (lA sup |a(t) -/3(i)| ) + e" T . (38) 

V 0<t<T J 

It easily follows from assumption (Ci) and from the dominated convergence theorem 
that 

limsupE[d(X x ,X x °)] < e~ T for every T > 0. 

Letting T -> +oo implies that lim x ^ Xo E[d(X x ,X x °)} = 0. 

(b) We deduce from (38) and from assumption (C2) that for every compact set K C R d , 
for every T > 0, 

limsup sup E[d(Y n ' x , X x )] < e~ T . 
Letting T -> +00 yields (36). □ 

Lemma 4. Assume that (r] n ) n >i and (7^) satisfy (C 3 e ) with e < 1 and (X£j. Then: 

(i) /or every t>0, for every bounded continuous function f : W d — > R, 

(n)/ r\ (n) I r\ n-*+oo n 

ft (w,/)-^ (w,/) — ► a.s.; 

(ii) ?/, moreover, a.s., every weak limit 1/ 00 ) (w, da) 0/ (i/^ n '(w,da)) n >i is i/ie dis- 
tribution of a Markov process with semigroup (Qt)t>o, then, a.s., i/ 00 -* (w, da) is £/ie 
distribution of a stationary process. 

Proof, (i) Let / :R d — > R be a bounded continuous function. Since x[ k ^ = Xr N , k t) , we 
have 

JH ' 1 fe=l 

From the very definition of N(n,T) and r(n,T), one checks that iV(fc — 1,T) < n — 1 if 
and only if r(n, T) > /c. Then, 

^ n ^ '"(n.t) 

TT^2 r lkf(Xr k _ 1 ) = — ^2 VN(k-i,t)+if(Xr N(k _ lt) ) 
" fe=i ™ fe=i 



1 " 

5-E 1J ^ Xr '-i) 1 { i - 1 S 



rr / v-- 1( ,-i/- l »-l^JV({0,...,n}i)}- 
™ fc=l 
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It follows that 



/{ t n) (u,f) - v { n) (oJ,f) = Tj~ Yl ( % ~ ? ?JV(fe-i,t)+i)/(^r w(fe _ 1 , t) ) 



fe=i 



1 " 

3- ^/(^(t-i.f)) 



r(n,t)+X 
1 " 

71 _ X] 7?fc ^^ r '=-i) 1 { fe - 1 ^ Ar ({ ---"}' t )}- 

™ fc=l 



Then, since / is bounded and since 



n n T"(">t) 

X] ? /fe 1 {fc-l^JV({0,. ..,«},{)} = _ Y2 ^(fc-l,i)+l 

fc=l fe=l fe=l 

r(n,t) n 

< X] l 7 ' fc ~ VN(k-l,t)+l\ + Vk, 
k=l k=T(n,t) + l 



we deduce that 



/ -, r{n,t) n 

\p\ n) (u>, /) - ( W) /)| < 2||/|U — ^ | % - VN{k . ht]+1 1 + — ^ 



ff n ^ "'l*-^;-!-*. £• 

\ fc=l fe=r(n,t)+l / 

Hence, we have to show that the sequences of the right-hand side of the preceding in- 
equality tend to 0. On the one hand, we observe that 

JV(fe-l,T) + l 2V(fc-l,T)+l 

\Vk ~ VN(k-i,t)+x\ < Y] |%-%-i|<max V j e . 

e=k+i - " £=fc 

Using the fact that X)fcfc _1T ^ +1 It — T + 71 and condition (14) yields 



r(n,t) 

JT l»7fc — »7JV(fc-i,*)+i| Tl ^— ?° 0. 



k=l 

On the other hand, by (C3 £ ), we have 



fc=r(n,T) + l " fe=r(n,r)+l 



k—r{n,T) + l 

which completes the proof of (i) 



Ik < TTl _ e — ► a.s. 
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(ii) Let Q + denote the set of non-negative rational numbers. Let (fe)e>i be an every- 
where dense sequence in Ck endowed with the topology of uniform convergence on 
compact sets. Since Q + and (fe)e>i are countable, we derive from (i) that there exists 
fl C such that P(f2) = 1 and such that for every u> G (l, every t G Q + and every I > 1, 

v\ (w,/e)-fo K.W — > 0. 
Let w G f2 and let (w, da) denote a weak limit of (i/' n '(w,da)) n >i. We have 

p[ oo) (u,h) = vi° o) (u,f l ) V*gQ+W>1 

and we easily deduce that 

^ (oo) (to, /) = i/<°°> ( W , /) Vi G M+ V/ G C A -(M d ). 

Hence, if z/°°) (u>, da) is the distribution of a Markov process (Y t ) with semigroup (Qf)t>o, 
we have, for all feC K (K. d ), 

J QtfiM^ (w, dx) = 1 /(x)^ («, da;) Vt > 0. 

Vq°°^ (u), dx) is then an invariant distribution for (Yt). This completes the proof. □ 



4.2. Proof of Theorem 1 

Thanks to Lemma 1(a) applied with £ = B(M + , M d ) and d defined by (6), 

i>V,da)^P, (da) a,.s.<=>v {n) (Lu,F) n ^° [ F(x)F„ (dx) a.s. (39) 



for every bounded Lipschitz functional F:B(R + ,R d ) — * K. Now, consider such a func- 
tional. By the assumptions of Theorem 1, we know that a.s., (vq (oj, dx))„>i converges 
weakly to v Q . Set $ F (x) := E[F(X X )], x G R d . By Lemma 3(a), $ F is a bounded contin- 
uous function on M. d . It then follows from (Co,i) that 



a.s. 



^-^^(X^Y^ [ $> F ( x )v Q (dx) = f F(x)W V0 (&x) 
n fc=1 J J 

Hence, the right-hand side of (39) holds for F as soon as 
i " 

±Y.Vk(F(X^y)-<P F (xt iy )) n ^0 a.s. (40) 



U k=l 
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Let us prove (40). First, let T > and let Ft be denned by (7). By Lemma 2, 

— J2 VkF T (X^) -— ]T Ws E[F T (X«<- 1 '>)/F rh _ 1 ] 



n — >+oo 



a.s. 



" fe=l " fc=l 

With the notation of Lemma 3(b), we derive from assumption (C2)(i) that 

E[F T (X^)/^ rk _ 1 } = ^ (4 k -V). 
Let N EN. On one hand, by Lemma 3(b), 



n 



n — >+oo „ 
y(*-l)|x«rl * U 8-S. 



fe=l 



On the other hand, the tightness of (v^ {lo, di')) n >i on W 1 yields 

1>(w, N) := sup(^" ) ( W> (5(0, iV) c ))) a.s. 



It follows that, a.s., 

supfi^^i^^r 13 )-^^- 15 )!^,^-) 



|>iV} 



<2||F|| oo V(« > JV) Ar =^ oo 0. 



Hence, a combination of (42) and (43) yields 



a.s. 



fc=i 



(41) 



(42) 



(43) 



(44) 



Finally, let (H*)^>i be a sequence of positive numbers such that, Ti — > +oo when £ — > +oo. 
Combining (44) and (41), we obtain that, a.s., for every l> 1, 



lim sup 

n — > + oo 



fc=l 



< lira sup 



1 " 

fin , 



fc=l 



+ lim sup 



n 



fc=i 
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By the definition of d, \F - F Te \ < e~ Tt ■ Then, a.s., 
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lim sup 

n — >+oo 



1 " 

1 J2v k (F(X^)-^(xt 1] )) 



H, 



k=l 



< 2c 



-T, 



W> 1. 



Letting I— > +oo implies (40). 

The generalization to non-bounded functionals in Theorem 1 is then derived from (28) 
and from a uniform integrability argument. 



4.3. Proof of Theorem 2 

(i) We want to prove that the conditions of Lemma 1(b) are fulfilled. Since (v^ (oj, dx))„>i 
is supposed to be a.s. tight, one can check that for every bounded Lipschitz functional 
F:D(R + ,R d ) — >R, (40) is still valid. Then, let {Ft)i>\ be a sequence of bounded Lipschitz 
functionals. There exists tl C ft with F(fl) = 1 such that for every lu G f2, (y^ (u>, dx)) n >i 
is tight and 

^^Vk(Fi(X^- 1 \u))-9 Ft {X^ k - 1 \u))) r ^ro W>1. (45) 

k— 1 

Let lu e & and let 4> u :N^ N be an increasing function. As (z/q ' n ^'(a;, dx))„>i is tight, 
there exists a convergent subsequence (v^""^"^^, dx)) n >i- We denote its weak limit 
by Voo. Since $ Fe is continuous for every £ > 1 (see Lemma 3(a)), 

We then derive from (45) that for every I > 1 

i/W«°*«( n )>(w,fi) ' l=± >°° /" ff(a)P^ (da). 
It follows that the conditions of Lemma 1(b) are fulfilled with U = {P^, /i G X}, where 
1= \ fj,eV(R d ),3cueQ and an increasing function 0:N^N,/i= lim v^ n)) (lu, da) I. 

n — > + oo 

Hence, by Lemma 1(b), we deduce that (v^ n > (lu, da))„>i is a.s. tight with W-valued limits. 
Finally, Theorem 2(ii) is a consequence of condition (14) and Lemma 4(h). 
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5. Path-dependent option pricing in stationary 
stochastic volatility models 

In this section, we propose a simple and efficient method to price options in stationary 
stochastic volatility (SSV) models. In most stochastic volatility (SV) models, the volatil- 
ity is a mean reverting process. These processes are generally ergodic with a unique 
invariant distribution (the Heston model or the BNS model for instance (see below) but 
also the SABR model (see Hagan et al. [8]), . . .). However, they are usually considered 
in SV models under a non-stationary regime, starting from a deterministic value (which 
usually turns out to be the mean of their invariant distribution). However, the instanta- 
neous volatility is not easy to observe on the market since it is not a traded asset. Hence, 
it seems to be more natural to assume that it evolves under its stationary regime than 
to give it a deterministic value at time 0. 4 

From a purely calibration viewpoint, considering an SV model in its SSV regime will 
not modify the set of parameters used to generate the implied volatility surface, although 
it will modify its shape, mainly for short maturities. This effect can in fact be an asset 
of the SSV approach since it may correct some observed drawbacks of some models (see, 
e.g., the Heston model below). 

From a numerical point of view, considering SSV models is no longer an obstacle, es- 
pecially when considering multi-asset models (in the unidimensional case, the stationary 
distribution can be made more or less explicit like in the Heston model; see below) since 
our algorithm is precisely devised to compute by simulation some expectations of func- 
tionals of processes under their stationary regime, even if this stationary regime cannot 
be directly simulated. 

As a first illustration (and a benchmark) of the method, we will describe in detail 
the algorithm for the pricing of Asian options in a Heston model. We will then show 
in our numerical results to what extent it differs, in terms of smile and skew, from the 
usual SV Heston model for short maturities. Finally, we will complete this section with 
a numerical test on Asian options in the BNS model where the volatility is driven by a 
tempered stable subordinator. Let us also mention that this method can be applied to 
other fields of finance like interest rates, and commodities and energy derivatives where 
mean-reverting processes play an important role. 



4 Whcn one has sufficiently close observations of the stock price, it is in fact possible to derive a rough 
idea of the size of the volatility from the variations of the stock price (see, e.g., Jacod [10]). Then, using 
this information, a good compromise between a deterministic initial value and the stationary case may 
be to assume that the distribution (iq of the volatility at time is concentrated around the estimated 
value (see Section 2.2 for application of our algorithm in this case). 
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We consider a Heston stochastic volatility model. The dynamic of the asset price process 
(St)t>o is given by 5*0 = so and 

dS t = S t (rdt + y/(l-pZ)vtdW? + pyfcdW?), 
dv t =k(9-v t )dt + ^ t dW 2 , 

where r denotes the interest rate, (W^ 1 ,^ 2 ) is a standard two-dimensional Brownian 
motion, p G [—1,1] and k, 9 and <; are some non-negative numbers. This model was 
introduced by Heston in 1993 (see Heston [9]). The equation for (vt) has a unique (strong) 
pathwisc continuous solution living in R + . If, moreover, 2k9 > then (vt) is a positive 
process (see Lamberton and Lapeyre [11]). In this case, (vt) has a unique invariant 
probability Vq. Moreover, Vo = 7(0, b) with a = (2k) /^ 2 and b = (2k9)/s 2 . In the following, 
we will assume that (vt) is in its stationary regime, that is, that 



5.1.1. Option price and stationary processes 

Using our procedure to price options in this model naturally needs to express the option 
price as the expectation of a functional of a stationary stochastic process. 

Naive method, (may work) Since (vt)t>o is stationary, the first idea is to express the 
option price as the expectation of a functional of (vt)t>o- by Ito calculus, we have 



St = so exp Mrt-iJ Vs ds} + p J ^T s dW 2 + sfl-p 2 J V^dW} 



(46) 



_ vt — vo — k6t + k J Q v s ds 



Since 

/ V£dW? = A(t,(«t)) 
Jo s 

it follows by setting M t = J ^Jv^dWl that 

S t = *(t,(v s ),(M s )), (47) 
where ^> is given for every t > 0, u and w € C(R+,R) by 

&(t,u,w) = s expert - \ J u(s) ds^j + pA(t, u) + ^/l- p 2 w(t)^j . 



Then, let F:C(R+,R) — > R be a non-negative measurable functional. Conditioning by 
yields 

M[F T ((S t ) t > )}=E[F T ((v t ) t > )}, 
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where, for every u G 

F T (u)=E 



F T ( (* (t, u, jf dW?) ) ) 



For some particular options such as the European call or put (thanks to the Black- 
Scholes formula), the functional F is explicit. In those cases, this method seems to be 
very efficient (see Panloup [20] for numerical results). However, in the general case, the 
computation of F will need some Monte Carlo methods at each step. This approach is 
then very time-consuming in general - that is why we are going to introduce another 
representation of the option as a functional of a stationary process. 

General method, (always works) We express the option premium as the expectation 
of a functional of a two-dimensional stationary stochastic process. This method is based 
on the following idea. Even though (vt,Mt) is not stationary, (St) can be expressed as a 
functional of a stationary process (vt,yt)- Indeed, consider the following SDE given by 



dy t = -vtdt + y/vidW}, 

dv t = k(6-v t )dt + <;jv- t dW?. 



(48) 



First, one checks that the SDE has a unique strong solution and that assumption (Si) is 
fulfilled with V(x±,X2) = 1 + xf -\-x\. This ensures the existence of an invariant distribu- 
tion vq for the SDE (see, e.g., Pages [17]). Then, since (vt) is positive and has a unique 
invariant distribution, the uniqueness of the invariant distribution follows. Then, assume 
that C(y 0l v Q ) = 9 . Since (v t ,M t ) = (v t ,y t - y + J Q y s ds), we have, for every positive 
measurable functional F : 



nF T ((S t ) t > )]=E[F T (Mt,v t ,M t )) t > )] 



F T [ [i>[t,vt,y t -yo+ / y s ds 

/ / t>0 



(49) 



where Py is the stationary distribution of the process (vt,yt)- Every option price can 
then be expressed as the expectation of an explicit functional of a stationary process. We 
will develop this second general approach in the numerical tests below. 



Remark 9. The idea of the second method holds for every stochastic volatility model 
for which (St) can be written as follows: 




where, for every i e {1, . . . hi :R+ — > M is a positive function such that hi(x) = o(|x|) 
as — ► +00, (Y t l ) is a square-integrable centered Levy process and (vt) is a mean 
reverting stochastic process solution to a Levy driven SDE. 
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In some complex models, showing the uniqueness of the invariant distribution may be 
difficult. In fact, it is important to note at this stage that the uniqueness of the invariant 
distribution for the couple (vt,yt) is not required. Indeed, by construction, the local 
martingale (M t ) does not depend on the choice of yo. It follows that if £(yo,vo) = fi, 
with fj, constructed such that C(v ) = vq, (49) still holds. This implies that it is only 
necessary that uniqueness holds for the invariant distribution of the stochastic volatility 
process. 

5.1.2. Numerical tests on Asian options 

We recall that (v t ) is a Cox-Ingersoll-Ross process. For this type of processes, it is well 
known that the genuine Euler scheme cannot be implemented since it does not preserve 
the non- negativity of the (vt). That is why some specific discretization schemes have 
been studied by several authors ( Alfonsi [1] , Deelstra and Delbaen [5] and Berkaoui et al. 
[4, 6]). In this paper, we consider the scheme studied by the last authors in a decreasing 
step framework. We denote it by (vi). We set Vq = x > and 

wr„ +1 =|ur n +fc7n+i(6 , -ur„)+cV^rT(^r„+i-^r„)l- 
We also introduce the stepwise constant Euler scheme (y t ) of (yt)t>o defined by 

yr„ +1 = i7r„ - 7n+iJ/r„ + y/^(Wr n+1 ~ ^rj. m = V G K d - 

Denote by {v[ ) and (fit) the shifted processes defined by Uj := vr k +t and y\ k ^ = 
yr k +t, and let (y' n '(w,da)) n >i be the sequence of empirical measures defined by 

1 ™ 

v {n) {uj,Aa) = — ^%l {(s ( fc -i) iS ( fc -i ))edQ} . 

n k=l 

The specificity of both the model and the Euler scheme implies that Theorems 1 and 2 
cannot be directly applied here. However, a specific study using the fact that (9) holds 
for every compact set of W + x K when 2k6/q 2 > 1 + 2v6/? (see Theorem 2.2 of Berkaoui 
et al. [4] and Remark 9) shows that 

v {n) (u J ,da) n ^ 00 Po (da) a.s. 

when 2k6/^ 2 > 1 + 2V6/?. Details are left to the reader. 

Let us now state our numerical results obtained for the pricing of Asian options with 
this discretization. We denote by C as (vQ,K,T) and P as {vQ,K,T) the Asian call and put 
prices in the SSV Heston model. We have 



C as ( I / ,^,T)=e- rT E 



mi 



1 f 
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and 



P as (v ,K,T) = c- rT E, 



K-- I 



With the notation of (49), approximating C as (vQ,K,T) and P as (yo,K,T) by our proce- 
dure needs to simulate the sequences (C™ s )n>i and (P™ s ) n >i defined by 



a 



i n 

m _ J_ \ " 
as tt Z^/ 



,-r-T 



k=l 



1 

n fc=i 



,-r-T 



- / *(s,w ( ' £ - 1) ,M( fe - 1 ))d S - A' 
^ Jo 



These sequences can be computed by the method developed in Section 1.3. Note that 
the specific properties of the exponential function and the linearity of the integral imply 
that (J *(i,w( n - 1 ',M("- 1 ))ds) can be computed quasi-recursively. 
Let us state our numerical results for the Asian call with parameters 



s = 50, 
= 0.01, 



r = 0.05, 

<r = 0.1, 



T=l, 
k = 2. 



p = 0.5, 



(51) 



We also assume that K £ {44, . . . , 56} and choose the following steps and weights: 7„ = 
rj n = n -1 / 3 . In Table 1, we first state the reference value for the Asian call price obtained 
for N = 10 s iterations. In the two following lines, we state our results for N = 5.10 4 and 
N = 5.10 5 iterations. Then, in the last lines, we present the numerical results obtained 



Table 1. Approximation of the Asian call price 



K 


44 


45 


46 




47 


48 


49 


50 


Asian call (ref.) 


6.92 


5.97 


5.04 




4.12 


3.25 


2.46 


1.78 


N = 5 ■ 10 4 


6.89 


6.07 


5.07 




4.13 


3.18 


2.49 


1.77 


N = 5 ■ 10 5 


6.90 


6.02 


5.00 




4.11 


3.24 


2.46 


1.79 


N = 5 ■ 10 4 (CP parity) 


6.92 


5.96 


5.04 




4.13 


3.26 


2.46 


1.78 


N = 5 ■ 10 5 (CP parity) 


6.92 


5.97 


5.04 




4.12 


3.25 


2.47 


1.78 


K 


51 


52 




53 




54 


55 


56 


Asian call (ref.) 


1.23 


0.82 




0.53 




0.33 


0.21 


0.12 


N = 5 ■ 10 4 


1.21 


0.81 




0.51 




0.34 


0.22 


0.11 


N = 5 ■ 10 5 


1.23 


0.82 




0.53 




0.33 


0.21 


0.13 


N = 5 ■ 10 4 (CP parity) 


1.23 


0.82 




0.53 




0.31 


0.21 


0.12 


AT = 5-10 5 (CP parity) 


1.23 


0.82 




0.53 




0.33 


0.21 


0.13 
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using the call-put parity 

C as {u , K, T) - P as (v ,S Q , K, T) = - e" rT ) - K^ rT (52) 

rl 

as a means of variance reduction. The computation times for N = 5.10 and N = 5.10 5 
(using MATLAB with a Xeon 2.4 GHz processor) are about 5 s and 51 s, respectively. In 
particular, the complexity is quasi-linear and the additional computations needed when 
we use the call-put parity are negligible. 

5.2. Implied volatility surfaces of Heston SSV and SV models 

Given a particular pricing model (with initial value so and interest rate r) and its asso- 
ciated European call prices denoted by C enr (K, T), we recall that the implied volatility 
surface is the graph of the function (K, T) i— ► a lmv {K, T), where cri mp (A', T) is defined for 
every maturity T > and strike K as the unique solution of 

C BS (so,K, T, r, a imp (K, T)) = C euT (K,T), 

where CBs(so,K,T,r,a) is the price of the European call in the Black-Scholes model 
with parameters So, r and a. When C eul (K,T) is known, the value of a- lmp (K, T) can be 
numerically computed using the Newton method or by dichotomy if the first method is 
not convergent. 

In this last part, we compare the implied volatility surfaces induced by the SSV and SV 
Heston models where we suppose that the initial value of (vt) in the SV Heston model is 
the mean of the invariant distribution, that is, we suppose that vq = 0. 5 We also assume 
that the parameters are those of (51), except the correlation coefficient p. 

In Figures 1 and 2, the volatility curves obtained when T = 1 arc depicted, whereas in 
Figures 3 and 4, we set the strike K at K = 50 and let the time vary. These representations 
show that when the maturity is long, the differences between the SSV and SV Heston 
models vanish. This is a consequence of the convergence of the stochastic volatility to its 
stationary regime when T — * +oo . 

The main differences between these models then appear for short maturities. That is 
why we complete this part by a representation of the volatility curve when T = 0.1 for 
p = and p = 0.5 in Figures 5 and 6, respectively. We observe that for short maturities, 
the volatility smile is more curved and the skew is steeper. These phenomena seem 
interesting for calibration since one well-known drawback of the standard Heston model 
is that it can have overly flat volatility curves for short maturities. 

5.3. Numerical tests on Asian options in the BNS SSV model 

The BNS model introduced in Barndorff-Nielsen and Shcphard [3] is a stochastic volatility 
model where the volatility process is a Levy-driven positive Ornstein-Uhlcnbcck process. 



5 This choice is the most usual in practice. 
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Figure 1. p = 0, K i-t- a ijnp (K, 1). 



The dynamic of the asset price (St) is given by St = So exp(X t ), 



dX t = (r- ±v t )dt + ^r t dWt+pdZ t , p<0, 
dv t = — fiv t dt + dZ t , [i > 0, 




Figure 2. p = 0.5, <r imp (lf, 1). 
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nuriufily 



Figure 3. p = 0, T ^ a imp (50,T). 



where {Z±) is a subordinator without drift term and Levy measure tt. In the following, 
we assume that (Z t ) is a tempered stable subordinator, that is, that 



7r(dy) = l {a>0 } 



cexp(-Ay) 



y 



l+a 



dy, c>0,A>0,ae (0, 1). 



As in the Heston model, we want to use our algorithm as a way of option pricing when 
the stochastic volatility evolves under its stationary regime and test it on Asian options 
using the method described in detail in Section 5.1. This model does not require a specific 



f 0.096 




-tfi — fir — 4r— o — a- ■ & ■ ta ■ 



12 3 



Figure 4. p = 0.5, T f— » crimp (50, T) . 
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discretization and the approximate Euler scheme (P) (see Section 3.2) relative to (vt) 
can be implemented using the rejection method. In Table 2, we present our numerical 
results obtained for the following choices of parameters, steps and weights: 

p=-l, A = ju = 1, c = 0.01, a=|, In =Vn = n~ 1/3 . 

The computation times for N = 5.10 4 and N = 5.10 5 are about 8.5 s and 93 s, respectively. 
Note that for this model, the convergence seems to be slower because of the approximation 
of the jump component. 

0.16 
S is 

on 

0.13 
t 6 0.12 
0.11 
0.1 



MM 

Figure 6. p = 0.5, Ti-f <7 imp (50,T). 
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K 44 45 46 47 48 49 50 



Asian call (ref.) 6.75 5.83 4.93 4.05 3.18 2.35 1.57 

Ar = 5-10 4 6.83 5.91 5.01 4.10 3.22 2.35 1.51 



N = 5 ■ 10 5 

N = 5 ■ 10 4 (CP parity) 
N = 5 ■ 10 5 (CP parity) 


6.78 
6.76 
6.75 


5.86 
5.85 
5.83 


4.96 
4.94 
4.93 




4.06 
4.07 
4.04 


3.19 
3.20 
3.17 


2.34 
2.29 
2.32 


1.52 
1.51 
1.54 


K 


51 


52 




53 




54 


55 


56 


Asian call (ref.) 
N = 5 ■ 10 4 
N = 5 ■ 10 5 

N = 5 ■ 10 4 (CP parity) 
Af = 5-10 5 (CP parity) 


0.91 

0.77 
0.79 
0.79 
0.83 


0.55 

0.46 
0.48 
0.47 
0.50 




0.39 

0.33 
0.34 
0.37 
0.36 




0.29 

0.27 
0.27 
0.27 
0.28 


0.23 

0.22 
0.21 
0.23 
0.22 


0.18 

0.19 
0.17 
0.19 
0.17 
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